module Equations
use sysutils
use strutils


contains


function SafeLn(Rate)
   implicit none
   double precision :: SafeLn
   double precision, intent(in) :: Rate
   if (Rate <= 0) then
      SafeLn = 0
      !call MessageOut("Error: Zero or negative argument to Log")
      return
   end if
   SafeLn = Log(Rate)
end function SafeLn


function SafeExp(Arg)
   implicit none
   double precision :: SafeExp
   double precision, intent(in) :: Arg
   if (87 < Arg) then
      SafeExp = 0
      !call MessageOut("Error: Argument to Exp > 87")
      return
   end if
   SafeExp = Exp(Arg)
end function SafeExp


function StirlingApprox(z)
!   Purpose: Ln to Gamma approximation
!   In:      ObsVal - observed value
!            ModVal - modified value
!   History: When        Who        What
!            15-08-07    Skalevik   Created
   implicit none
   double precision, intent(in) :: z
   double precision :: StirlingApprox
   
   if(z > 0) then
      StirlingApprox = 0.5 * (1.8378771 - log(z)) + z * (log(z + 1/(12.0 * z - 1 / (10.0 * z))) - 1)
   else
      StirlingApprox = 0
   endif
end function StirlingApprox
 

function rndn(rintvar)
! returns a random variable between 0 and 1, adjusts the
! seed rintvar
   double precision, intent(inout) :: rintvar
   double precision :: rndn
   rintvar = rintvar * 23.0
   rintvar = mod(rintvar, dble(100000001.0))
   rndn = rintvar / 100000001.0
end function rndn


function snrn(r)
! WHen r is a random number between 0 and 1
! the function returns the number with that probability
! in a standard normal distribution
   implicit none
   double precision snrn, r, t, sg
   if (r > 0.5d0) then
      t = sqrt(log((1.0d0 / (r - 0.5d0))**2))
      sg = 1.0d0
   else if (r < 0.5) then
      t = sqrt(log((1.0d0 / (r - 0.5d0))**2))
      sg = -1.0d0
   else
      sg = 0.0d0
      t = 0.0d0
   endif
  ! Smith and Bain approximation to std normal distribution:
   snrn = sg * (t - (2.30753d0 + 0.27061d0 * t) / (1.0d0 + 0.99229d0 * t + 0.04481d0 * t * t))
end function snrn

      
end module Equations
